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Abstract 

We study an improved three-dimensional Ising model with external magnetic field 
near the critical point by Monte Carlo simulations. From our data we determine 
numerically the universal scaling functions of the magnetization, that is the equation 
of state, of the susceptibility and of the correlation length. In order to normalize 
the scaling functions we calculate the critical amplitudes of the three observables 
on the critical line, the phase boundary and the critical isochore. The amplitudes 
lead to the universal ratios C+/C- = 4.756(28), = 1.723(13), Qc = 0.326(3) and 
Q2 = 1.201(10). We find excellent agreement of the data with the parametric repre- 
sentation of the asymptotic equation of state as found by field theory methods. The 
comparison of the susceptibility data to the corresponding scaling function shows a 
marginal difference in the symmetric phase, which can be explained by the slightly 
different value for used in the parametrization. The shape of the correlation- 
length-scaling function is similar to the one of the susceptibility, as expected from 
earlier parametrizations. The peak positions of the two scaling functions are coin- 
ciding within the error bars. 
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1 Introduction 



The Ising model is the simplest non-trivial spin model with critical behaviour. In 
three dimensions the universality class of the Ising model includes as well many real 
physical systems as also theoretical models which exhibit a continuous phase transi- 
tion. Important examples are known from condensed matter physics and high-energy 
physics (see for instance the detailed review in Ref. [T]). Due to their simplicity and 
the many physically relevant applications the 3d Ising model and members of the 
corresponding universality class belong to the most extensively studied systems. 
In these investigations both a variety of analytic methods P]-|S1 and Monte Carlo 
simulations have been used to calculate the universal quantities characterizing the 
universality class, such as critical exponents and amplitude ratios. 

In the vicinity of the critical point the behaviour of the order parameter, the mag- 
netization, and the correlation length are described by asymptotic scaling functions. 
The functions depend on the temperature and on the magnetic field and are univer- 
sal apart from metric factors which are specific for the model under consideration. 
They are asymptotic in so far as they are valid in the thermodynamic limit, that 
is for V —>■ oo, and for T — * T^, so that irrelevant scaling fields can be neglected. 
In order to be comparable to other model systems the scaling functions have to 
be normalized which amounts to fixing the metric factors. There exist numerous 
parametrizations (see e. g. Refs. ^ and |lI-[Zj) of the equation of state, the scaling 
function of the order parameter, and other scaling functions. Since a parametric 
form of a scaling function has to describe the critical behaviour of the respective 
observable in accord with the known critical exponents and amplitude ratios these 
are usually taken either as an input to fix the parameters or as a check. Further 
critical quantities may then eventually be calculated from the parametrizations. 

Instead of relying only on the known critical parameters one can in principle 
directly determine the scaling functions from simulation data, taking at most the 
critical exponents as input. This has already been achieved in the case of the three- 
dimensional 0(2) |8| and 0(4) vector models 0110] for the equation of state. In the 
3d Ising model class (the 0(1) vector model) no such attempt has been undertaken 
as far as we know. In this paper our main intention is to accomplish the same for 
that universality class and to compare the parametric representation of the equation 
of state as found from field theory directly to simulation results. A further aim 

is to find numerically the scaling function of the correlation length. It is clear, that 
for this project primarily simulations with non-zero magnetic field are required. The 
usual simulations at zero magnetic field can at best serve to normalize the scaling 
functions. The principal difficulties which we expect are due to the use of finite 
volumes in the simulations and to the finiteness of the critical region where scaling 
works. In order to attenuate these problems we use relatively large lattices with up 
to 120 points in each direction and we employ an improved Ising model to reduce 
corrections to scaling. 
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The model which we want to study is the one-component 0^ or Landau- Ginzburg 
modeL It is defined by 

= -J 0.0, + YM>l + K<t>l-if] - hY,^,- (1) 

<x,y> X X 

Here x and y are the nearest-neighbour sites on a three-dimensional simple cubic 
lattice, (px is the field variable at site x and H is the external magnetic field. We 
consider the coupling constant J as inverse temperature, that is J = 1/T. In contrast 
to the ordinary Ising model the length of the spins (px in Eq- © is not fixed to one. 
For A = one obtains the Gaussian model, whereas in the limit \ ^ oo the simple 
Ising model is recovered. By choosing an appropriate A-value it is then possible to 
eliminate leading order corrections to scaling. This has been shown numerically in 
Refs. jTT]-P3] for H = 0. There is however a drawback in using a finite A: we cannot 
directly verify our results by comparison with the known, non-universal amplitudes 
(see e. g. Ref. of the normal Ising model. 



The observables which we want to measure in our simulations are the magneti- 
zation M, the susceptibility x &nd the correlation length ^. The magnetization is 
the expectation value of the lattice average of the spin variable 

M = il^Y^x) = (0), (2) 

X 

where V = L"^ and L is the number of lattice points per direction. As long as H 
is non-zero and the volume large enough there is no problem with this definition. 
At zero magnetic field, H = 0, however, the lattice average of the spins will have 
a vanishing expectation value on all finite lattices, that is ( ) = . At very small 
magnetic fields and/or small L we approximate M therefore by J3] 

M ^ (101) . (3) 

On finite lattices the magnetization of Eq. Q approaches the infinite volume limit 
from above, whereas M as defined by Eq. Q for if 7^ reaches the thermodynamic 
limit from below. The susceptibility is defined as usual by the derivative of the 
magnetization 

X = ^ = V{{<P')-M^). (4) 

For the determination of the correlation length we consider three different definitions. 
The second moment correlation length is obtained from 

^ 1/2 

f^2 ' 



where the Hn are the moments of the connected part of the correlation function 

Hn = ^a;"(0.0o)c, (6) 

X 

(0x00 )c = (0x0o)-(0)', (7) 



2 



and a;^ = x^. The moment iiq is equal to the susceptibihty. A more and more often 
used definition of the correlation length is 



where F is the Fourier transform of the correlation function at momentum — 
27re^/ L, and a unit vector in one of the three directions 

^ = ^{\J2^Mv)<p.\') ■ (9) 

X 

In the simulations we compute F from an average over all three directions. Finally, 
the exponential or "true" correlation length ^^xp governs the large distance behaviour 
of the connected part of the correlation function 

( 0x00 )c ~ exp(-|x|/^exp) ■ (10) 

We explain later in detail how we determine this observable from the correlation 
function data. 

The rest of the paper is organized as follows. First we discuss the critical be- 
haviour of the observables and the universal scaling functions, which we want to 
calculate. In Section 3 we describe some details of our simulations. Section 4 serves 
to determine the amplitudes of the magnetization and the susceptibility. In the fol- 
lowing Section 5 we compare our data for the equation of state to a parametrization 
from field theory. The determination of the correlation length and its amplitudes is 
described in Section 6. In Section 7 we present our results for the scaling function 
of the correlation length. We close with a summary and the conclusions. 



2 Critical Behaviour and Scaling Functions 

In the thermodynamic limit {V oo) the observables show power law behaviour 
close to Tf.. It is described by critical amplitudes and exponents of the reduced 
temperature t = {T — T^jTc- The scaling laws at if = are for: 

the magnetization 

M = B{-tf for i < , (11) 

the susceptibility 

X = C^\t\-^ for t^±0 , (12) 
and the correlation length (without fixing for the moment a special definition) 

^ = for i ^ ±0 . (13) 
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On the critical line T = or t = we have for if > the scaling laws 

M = B^H^'^ oi H = D^M^ , (14) 

X = C^H^I^-^ with = , (15) 

and for the correlation length 

e = i'H-^^ , V, = u/pS . (16) 

We assume the following hyperscaling relations among the critical exponents to be 
valid 

2-a = di^, 7 = /3(5-l), di^ = /3(l + 5). (17) 

As a consequence only two critical exponents are independent. Because of the hyper- 
scaling relations and the already implicitly assumed equality of the critical exponents 
above and below one can construct a multitude of universal amplitude ratios [TE] 
(see also the discussion in Ref. PP). The following list of ratios contains those which 
we will determine here 

U2 = C+/C- , = t/r , (18) 

= C+D,B'~^ , = B\ey/C^ , (19) 

and 

Q2 = ie/ty^''c+/c' . (20) 

The critical behaviours of the magnetization and the susceptibility originate from 
the singular part of the free energy density. In the thermodynamic limit it obeys 
the scaling law 

fsiut,Uh) = h-'^fsihy^UuW^Uh) . (21) 
Here 6 is a free length rescaling factor, and 

yt = l/^ , Vh = l/^c . (22) 

We have neglected possible dependencies on irrelevant scaling fields Ui with expo- 
nents Ui < 0. Close to Tc and for if — »• the remaining relevant scaling fields Ut 
and Uh are proportional to the reduced temperature, Ut ~ t, and the magnetic field, 
Uh ~ if, that is one can replace Ut and Uh by a normalized reduced temperature 
t = ITc/Tq and magnetic field h = H/Hq. To be definite we assume in the following 
that if is positive. Choosing the scale factor b such that W^^h = 1 we find 

/, = h'^-^^.a/h^"'') , (23) 

and using M = —dfs/dH we arrive at 

M = h^'^fcit/h^'^^) , (24) 
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where and fc are universal scaling functions once the normalizations are fixed. 
The relation (j^^ is one form of the magnetic equation of state. Alternatively, one 
may use the Widom- Griffiths form of the equation of state ^Tj QBl 

y = fix) , (25) 

where 

y = h/M\ x = i/M^'^ . (26) 
We take the standard normalization conditions 

/(O) = 1 , /(-I) = , (27) 

which imply 

M{t = 0) = h^'^ and Hq = , (28) 

M{h = 0) = {-i f and To = E^^/^Te . (29) 

The two forms ()24|) and are of course equivalent. The function fdz) and its 
argument z are related to x and y by 

= y-V5 ^ ^ = t/h^/^^ = xy-^l^^ . (30) 

Correspondingly the normalization conditions ^Tj\ translate into 

/g(0) = 1 , and fa{z);^^J,-zf . (31) 



Since the susceptibility x is the derivative of M with respect to H we obtain from 
Eq. (El 

dM _ h^/^-^ 

'oh ~ ~irn 

with 

I f „ , . z 



X = TTTF = ^^/x(^) . (32) 



/x(^) = ^ (^/c(^) - ^/g(^) J • (33) 

The asymptotic behaviour of for if at fixed t, that is for z ±oo, is 
determined by Eq. (fT^ 

/x(^)_loo^^^^^'"'(±^)"'- (34) 

We note, that the amplitude for z oo is simply the universal ratio i?^, whereas 
for z — >• — oo it is the analogous quantity R~ = R^/U2. The leading terms for fc 
are respectively 

fciz) {-zr + R-i-z)-\ (35) 

fci^) ^loo ^x^""' (36) 



in accord with the normalization ()3H1 . The asymptotical identity fc = fx for -2 ^ C)0 
is due to the fact that for T > and small magnetic field M is proportional to H. 
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The correlation length fulfills a scaling law analogous to the one of the free energy 
density 

= bab'^^Ut,by'^UH) . (37) 

In the same way as before we obtain from this equation the dependence of the 
correlation length on t and h in the critical region and the thermodynamic limit 

e = h-''^g^{z) , (38) 

in terms of a scaling function g^{z). Here, we have assumed that the metric factors 
in the scaling functions for fg and ^ are the same, which is the usual assumption 
[T6] . The function g^{z) is then universal except for one normalization factor. On 
the critical line t = or z = we find from (fTH|l 

g^{0) = CD-/^ = e{Br'P , (39) 

and from (fT!?|) the asymptotic behaviour aX z ^ ±oo 

9d^) .=±oo^^B'^"^^^)"' ■ (40) 
Indeed, the ratio of the amplitude for 2; — > oo in (pn|) and (7^(0) is universal 



(41) 



whereas g^{0) itself is not universal. Pelissetto and Vicari 'J] associate a universal 
scaling function /g(x) to the correlation length by the relation 

^2 _ (5^)25/'iM-2-//5/5(x) , (42) 

where x is defined as in fl2fi|l . It is straightforward to show that 

The value x = corresponds to z = 0. Since both f^{x) and fci^) are universal this 
is also true for 

9d^) = 9d'^)lgm ■ (44) 

As we shall see, the scaling functions of the correlation length and of the suscep- 
tibility have similar shapes, each with a peak at a positive 2;-value. It is therefore 
useful to consider the ratio ^^/x- We can express it with the corresponding ratio of 
universal scaling functions 

^ = {^-f{B')-^H-'^^^ ■ . (45) 

Here rj = 2 — ^/1/ is the exponent describing the behaviour of the two-point function 
at the critical point. We postpone the discussion of parametrizations of the scaling 
functions to the sections where we compare them to our data. 
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3 Simulation Details 



In the numerical simulation of the model (0) we followed the examples given by 
Brower and Tamayo and Hasenbusch They combined a simple local 

Metropolis algorithm for the update of the field with a cluster algorithm, which 
only updates the sign of the fields (p^ at fixed values of the modulus \(f)x\- During the 
cluster update the system is therefore treated as an Ising model with link-dependent 
coupling constants J^y = J\(j)x\\(py\ and local magnetic fields = H\(j)x\. In order 
to achieve ergodicity of the entire update process the two algorithms have to be used 
in an alternating manner. In a local Metropolis step a new proposal 0^ for the field 
is produced by [13] 



where r is a uniformly in (0,1] distributed random number and s an adjustable 
size parameter. The prescription ()4(ij) obviously allows for changes of the lengths 
of the spins. For the cluster algorithm we chose Wolff's single cluster update ^0] . 
which was modified to include a non-zero magnetic field [21]. This variant was 
also successfully used in previous simulations of the 0(2) ^ and 0(4) spin models 
A complete Monte Carlo update consists then of n consecutive cluster updates 
followed by a Metropolis sweep of the entire lattice. 

As explained in the introduction we aim at minimizing the corrections to scal- 
ing. We have therefore fixed the parameter A of the model to the optimal value 
A = 1.1 found by Hasenbusch [12]. The very precise value Jc = 0.3750966(4) for 
the corresponding critical coupling was adopted from the same paper. All our sim- 
ulations were done on three-dimensional lattices with periodic boundary conditions 
and linear extensions ranging from L = 24 to L = 120. The size parameter s in 
Eq. fim|) was set to s = 2. That resulted in acceptance rates of 59% to 67% in the 
Metropolis step for the simulated J and if-values. The number n of cluster updates 
was chosen in the region 10 < n < 25 for T > Tc and in the cold phase (T < Tc) in 
the range 3 < n < 8. Between two measurements of the observables 20 to 50 Monte 
Carlo updates were performed, such that the integrated autocorrelation times of the 
energy and the magnetization were Tint ~ 3 for L < 48 and Tint ~ 5 for L > 72. In 
general we made 20000 measurements, for L = 24 sometimes 40000. 

The coupling constant region which we have explored was 0.365 < J < 0.4, the 
magnetic field was varied from H = to H = 0.005. The critical region, where 
asymptotic scaling works, is clearly inside this parameter range. We come back to 
this point when we discuss our data for the scaling functions. Further details on the 
Monte Carlo update, for example on the distribution of the modulus of the spins, 
can be found in Ref. 2^. 




(46) 
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4 The Critical Amplitudes of M and x 



We consider first our results for tlie magnetization at if = and T < Tc- Tlie 
primary objective liere is to determine tlie normalization Tq or equivalently the 
critical amplitude B of (jllll . In Fig. ^ we show M from Eq. (jH)) as a function of 
Tc — T for different lattice sizes. The temperatures correspond to coupling constants 
in the range 0.376 < J < 0.4 . In order to eliminate finite-size effects at a fixed J 
(or T), we have simulated at increasingly larger values of L. 

We have fitted the results from the largest lattices to the following ansatz taking 
into account possible analytic and non-analytic corrections to scaling 

M = 6o(Te - T)'3[l + 6i(T, - rr'^ + b2{T, - T)] . (47) 

For the fits we used the critical exponents of Ref. [3] as input 

P = 0.3258(14) , u = 0.6304(13) , u = 0.799(11) , (48) 

because we wanted to test the corresponding equation of state. The remaining ex- 
ponents were determined from the hyperscaling relations (fTTjl . Instead of these field 
theory results we could as well have used the exponents of Hasenbusch ^3] or the 
estimates proposed by Pelissetto and Vicari pp. The difference is marginal, the 
various exponent results coincide within the error bars. The fit with the ansatz ()47|) 
nicely confirms the expected absence of the leading non-analytic correction term: bi 
is zero within error bars. The best fit to the points with Tc — T < 0.10 and 6i = 
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Figure 1: The magnetization M at H = versus T^ — T for L = 24, 48, 72 and 96. 
The solid line shows the fit 
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Figure 2: The magnetization M at as a function of the magnetic field H for 
L = 24, 48, 72, 96 and 120. The sohd hue shows the fit (j53|. 



is given by 

M = 1.0735(2)(rc-T)^[l-0.061(2)(rc-r)] , 
with /Nf = 0.8 . From Bq we obtain 

To = 6o^/^ = 0.8044(4) , B = b^T^ = 1.4776(2) . 



(49) 



(50) 



On the critical isotherm, that is at Tc, we have measured the magnetization for 
> on lattices with L = 24,48,72,96 and 120 using Eq. 0. The results are 
shown in Fig. |21 We observe an increasing finite size dependence close to H = 0. 
The results from the largest lattices have been fitted to the ansatz 



M 



(51) 



In the i7-interval [0.00005, 0.0005], however, only the leading term is of relevance. 
A fit in this interval with simply the first term in ()5H) leads to 



B" = 1.0435(15) . 



(52) 



The remaining amplitudes BI2 niay then be determined from a subsequent fit to 
all available data. As on the coexistence line we find that the value for the leading 
non-analytic correction amplitude B^ is small and of the same size as its error. We 
have therefore made fits with fixed = in the region 0.00005 < H < 0.0025 and 
find 

M = 1.0435(15)/J^/^[1 - 2.65(5)//] . (53) 
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Figure 3: The susceptibility x at if = for L = 48, 72 and 96, (a) below Tc versus 
Tc — T and (b) above Tc versus T — Tc. The solid lines show the fits from Eqs. 
and (jnni), respectively. 



with x^/Nf ^ 2. From the value of B'^ in (|^^ we can immediately deduce the 
values of two further quantities 

= 0.2172(3) and = Hq = 0.8150(56) . (54) 

Our next aim is the determination of the critical amplitudes C^. To this end we 
have evaluated the data of the susceptibility below and above the critical tempera- 
ture for H = 0. The data points are plotted in Figs. El (a) and (b) as a function of 
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{Tf.—T) and {T—Tc), respectively. Inspired by our experience with the magnetization 
in the cold phase, we have directly used a fit ansatz without a non-analytic correc- 
tion term. The data points from the largest lattices in the range 0.01 < Tc — T< 0.10 
lead to 

X = 0.1006(5)(-t)-^[l - 2.0(3)(-t)] , (55) 

with /Nf ^ 1.3 . Note that here we have used the reduced temperature as variable. 
A corresponding ansatz in the symmetric phase (T > Tc) shows that the correction 
term is zero inside the error bars. The points from the largest lattices are however 
well described by the leading term only with 

X = 0.4785(15) , (56) 

Here, x^/^f ~ 1-2 . From the fits and (jSUj) we have 

C- = 0.1006(5) and C+ = 0.4785(15) . (57) 

This enables us to calculate the two universal amplitude ratios 

U2 = 4.756(28) and = 1.723(13) . (58) 

The ratio U2 is in excellent agreement with two other MC results: 4.75(3) from the 
standard Ising model |2S! and 4.72(11) from (3-|-l)-dimensional SU{2) gauge theory 
as well as other results (see Table 11 in P and Table 5 in jS]). The ratio is also 
well inside the bounds [4.688,4.891] derived in Ref. The result for is the 

first one solely from MC calculations. It is in agreement with 1.57(23) from Refs. 
[HIES! and slightly higher than the result 1.660(4) from Refs. fHITfj. 



5 The Magnetic Equation of State 



After having determined the normalizations Tq and Hq we can calculate the univer- 
sal scaling function of the magnetization for the Ising class. To this end we have 
performed further simulations at 14 couplings in the range 0.385 < J < 0.365 with 
external fields H varying between 0.0001 and 0.005 . We convinced ourselves that 
the used lattice extensions were large enough to produce volume independent re- 
sults. For that purpose L = 48 was sufficient at high temperatures (small J) , closer 
to Tc lattices with linear extension L = 72 and 96 were needed. The finite size ef- 
fects diminished with increasing field, generally they were stronger in the cold phase 
than in the hot phase. In Fig. |3]we present the results for faiz) obtained from the 
magnetization data in the coupling range 0.38 < J < 0.37. Obviously these data 
are scaling well. Since the data from the couplings J = 0.385 and 0.365 showed 
already slight deviations from scaling behaviour we have discarded them. Similarly, 
limitations were found for the scaling if-regions. The data for fciz) were scaling for 
J = 0.38 only up to = 0.0003, for J = 0.378 to J = 0.375 up to = 0.0005 and 
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Figure 4: The equation of state: fc = Mh^^/^ as a function of 2; = ih^^/^^ . The 
sohd hne shows the parametrization of Ref. [3 , the dashed hues the asymptotic 
forms and The numbers refer to the different J = 1/T- values of the data, 
the dotted hues cross at the normahzation point /g(0) = 1. 



for J = 0.3745 to J = 0.372 np to H = 0.0015. At J = 0.37 the scahng i7-range 
was hmited again hj H = 0.0003. Only data within these parameter ranges are 
shown in Fig. |3]. From the amplitude ratios U2 and in (I58j) one can calculate the 
asymptotic forms and of fci^)- As can be seen in Fig. IH the data for fc 
coincide with the corresponding asymptotic form already for z < — 4 and z>6~7. 

The numerically determined scaling function can now be compared to parametric 
representations of the equation of state. A prominent and widely accepted repre- 
sentation of the equation of state is that of Guida and Zinn- Justin |S|-|Sj, which is 
based on field theory methods. They parametrize the equation of state in terms of 
two new variables R and 6 by setting |28j-[Hn] 

M = moRf^e , (59) 
t = R{1- e^) , (60) 
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h = hoR^^h{e) , (61) 

where mo and ho are two normalization constants. Because M is linear in 6 this 
is also known as the linear model. The function h{6) is an odd function and was 
approximated by a 5th order polynomial with the following coefficients ^ 

h{e) = 6 - 0.76201(36)^3 + 0.00804(11)^^ . (62) 

In Refs. jUE] slightly different coefficients are given. The function h{6) in ()62|1 has 
two zeros. The one at 6' = corresponds to /i = 0, t > 0, the second at 9q = 1.154 
to the coexistence line h = 0, t < 0. For 6 = 1 one obtains the critical line t = 0. 
Using the normalization conditions (j27jl one finds 

^o= ^\ ' , ho = ^ (63) 
f^o h(l) 



The Widom- Griffiths scaling variables x and y are only functions of 9 

i-e^ (e„y" h(e) 



The same applies to the scaling function fa and its argument z 



= ^ T7TT • (65) 



In Fig. 0] the solid line represents the parametrization from Eq. (j65|l with h{9) from 
(jn2I)- Evidently, the agreement between the data and the solid line is excellent. 
Also, the asymptotic form for 2 < coincides directly with the parametrization in 
a large range, whereas the asymptotic form for 2; > is marginally higher in the 
shown large 2;-range. 

We have also plotted the data in the Widom- Griffiths form. The small x-region 
is shown in Fig. El (a), the large x-region in Fig. El (b). This kind of representation 
of the equation of state is asymmetric with respect to the critical point: the cold 
or broken phase (T < Tc) extends in x from —1 to and in y from to 1, the hot 
phase (T > T^) reaches in x from to infinity and in y from 1 to infinity. Since 
the magnetization appears here in both variables small deviations from scaling are 
easier to detect than in Fig.lH This applies in particular to the cold phase. Therefore 
we included only points for x < in Fig. ini(a) with H < 0.0003, because the points 
with larger if- values exhibited already very slight, but perceptible deviations. In the 
hot phase we could use all the points already displayed in Fig. jH As in Fig. (H the 
data show in the Widom-Griffiths representation a similarly impressive agreement 
with the parametrization of the equation of state by Guida and Zinn- Justin. 

The scaling function /^(-z) of the susceptibility is connected via Eq. (j33|) to the 
scaling function /g(^) of the magnetization. It is therefore in principal known. On 
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Figure 5: The Widom-Griffiths form of the equation of state: y = h/M^ as a function 
of X = i/M^/^ . The sohd hues show the parametrization of Ref. Pj, the numbers 
refer to the J = 1/T- values of the data. Part (a) of the figure displays the small 
x-region where the stars denote the normalizations (fTfj) . part (b) the large x-region. 



the other hand we have direct data for x from our simulations which allows for 
another check of parametric representations and the scaling hypothesis. In Fig. IHl 
we show the respective data for the same J- values as in Fig. |3J We observe explicit 
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Figure 6: The scaling function of the susceptibihty fxi^) = xHoh}~^/^ . The sohd 
hne is the parametrization of Ref. [Sj, the dashed hues are the asymptotic forms 
(jnH). The numbers refer to the J = 1/T- values of the data. The star represents the 
normalization fx{0) = 1/5- 

scaling of the data, sometimes even in a larger if-range. In comparison to the 
parametrization of Ref. we find complete agreement in the cold phase {z < 0) 
and a slightly increasing small difference to the data in the hot phase {z > 0). The 
reason for this little discrepancy lies in the different asymptotic amplitudes which are 
used in the parametrization and which we measured from our data. The amplitude 
in question is R^. In the parametrization scheme one derives 

Rx = hm/'{el~i)-^ . (66) 

This results in the number = 1.666 , whereas we found from our amplitude 
measurements on the critical line and at h = the value R^. = 1.723(13), that is a 
slightly higher number. The data shown in Fig. IHl are consistent with that finding: 
at larger 2;-values they coincide with the corresponding asymptotic curve. 



6 The Correlation Length 

Instead of using correlation functions of the individual spins it is more favourable 
|3T] to consider spin averages over planes and their respective correlation functions. 
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For example, the spin average over the (x, ?/)-plane at position z is defined by 

^^ = Jji^ • (67) 

The average S of all 5"^ is equal to the lattice average and 

(5,) = (5) = (0) . (68) 

Correspondingly, we define the plane-correlation function G(^z) by 

G(z)=L2((5o5,)-(5)2) . (69) 

Here, z is the distance between the two planes. Instead of choosing the 2;-direction 
as normal to the plane one can as well take the x- or ?/-directions. Accordingly, we 
enhance the accuracy of the correlation function data by averaging over all three 
directions and all possible translations. The correlators are symmetric and periodic 
functions of the distance r between the planes 

G{t) = G{-t) , and G{t) = G{L - r) . (70) 

The factor on the right-hand side of ensures the relation 

L-l L/2 

T=0 T=-L/2+l 

Similarly, the second moment fi2 may be expressed in terms of the plane correlators 
as 

L/2 

f,2 = d Yl ^'G{r). (72) 

r=-L/2+l 

Like the point-correlation functions in Eq. pO|) the plane correlators G{t) decay 
exponentially. In Fig.EJwe show the logarithm of typical plane-correlation functions 
for different L, here at T = Tc and H = 0.0003. The expected behaviour is actually 
confirmed, as well as the infiuence of periodic boundary conditions on finite lattices. 
In order to obtain the exponential correlation length from G{t) on a lattice with 
linear extension L and periodic boundary conditions we therefore make the ansatz 

G(r) = A[exp(-r/0+exp(-(L-r)/0] , (73) 

and then try to fit the data for G{t) in an appropriate r-range. The ansatz ()73|) 
implies of course, that there are no additional excitations contributing to G{t). In- 
spired by the experiences reported in [21] , we proceed in the following way. First we 
calculate an effective correlation length ^e//(T) from ()73|1 using only the correlators 
at r and r + 1. For r <^ L this correlation length is approximately given by 

^^"^''^ = ln(G(r+!)/G(r)) " ^^^^ 
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Figure 7: The plane-correlation function G{t) at T = T^. and H = 0.0003 for 
L = 48,72,96 and 120. The solid lines show the fits from the ansatz ()73|1 . the 
dashed lines indicate the r-range, inside which the fits were done. 



With increasing r also ^eff increases and eventually reaches a plateau inside its error 
bars. Most likely, the lower values of ^e// at small r are due to higher excitations. 
At large distances r the resulting ^e// start to fluctuate when the relative error of 
the data becomes too large. The two limits define an intermediate r-range where 
a global fit with the ansatz (fTHj) can finally be used to estimate the exponential 
correlation length. 

In order to give an impression of the actual implementation of this procedure we 
show in Tabled the results of the fits (^■^**) to the plane-correlation functions from 
Fig. The corresponding correlation length is here close to 8. The plateau of ^eff 
is reached for all L at a value for Tmin between 6 and 10. A reasonable upper fit 
limit is Tmax = 18 to 22. Varying Tmax in this range does not affect the fit result. 
Also a variation of from 6 to 10 leads only to changes inside the error bars. The 
respective x^/^f are less than or equal to 0.02. If we compare C,-^^^ to the effective 
correlation length at Tmin, we find essentially the same result, including the size of 
the error bar. 

We have calculated as well the correlation length C,^ from Eq. (jHl). We see 
from Table ^ that is compatible with but has a larger error. Only for the 
largest lattice with L = 120 we obtain a difference. This may be due to the use of all 
correlations in the calculation of the Fourier transform, also those at large distances, 
where the relative errors are large. In addition, the numerical determination of 

requires an increasing accuracy with increasing L because of the diminishing 
denomiator sin(7r/L) in Eq. (jH} and the difference x/F — 1 in the nominator. 
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Table 1: The correlation length at T = and H = 0.0003 from different methods: 
^fi.t from fits to the ansatz (f7S|) in the r-interval [Tmin,Tmax]] from Eq. (jH} and 
^2nd from Eq. (0) using Eq. ((72)). For comparison we show also (.effiTmin)- 

The accurate numerical calculation of ^2nd is also problematic. The reason for 
that is the weight in the sum of Eq. ()72j). which emphasizes the correlators 
G{t) with the largest relative errors. Caselle and Hasenbusch [2H] have therefore 
replaced the correlation function for r > Tmax ~ 3^ by its large distance form with 
^ = C,eff{Tmax)- The sum was then extended to infinity. The procedure provides 
an estimate of the infinite volume value of fi2- On the other hand the finite size 
dependence of C,2nd is then lost and the determination of ^2nd is no longer indepen- 
dent of the other definitions. Because of that, we have derived the second moment 
correlation length directly from Eq. ()72|1 . truncating the sum when the numerical 
correlation function result was compatibel with zero. The corresponding values of 
C,2nd in Tabled have relatively large errors. The values increase with increasing L, 
approaching the results from other definitions from below. At L = 120 the estimate 
for ^2nd is probably too low, because of the mentioned numerical difficulties. 

Our test of the different definitions of ^ was performed at T = Tc . Already here 
we have found indications for the presence of higher states. In general, the infiuence 
of higher exitations becomes even stronger in the broken phase. An extensive dis- 
cussion of this point and the actual determination of the first two higher states with 
a variational technique can be found in Ref . |S2] • In the same paper and in Ref . [21] 
the numerical estimates of the exponential correlation length in the broken phase 
have been compared to high precision results from the Z2 gauge model jHS], which 
is related to the Ising model by a duality transformation. The corresponding values 
differed at most in the third digit and showed that ^g// slightly underestimates the 
exponential correlation length. In the high temperature phase no nearby exitations 
are present and C,eff is a good estimator of the true correlation length. As a result of 
these considerations we resort in the following to the estimation of the exponential 
correlation length via the fit method. The corresponding scaling function can be 
calculated to a good approximation in this way. The extreme precision aimed at 
in Refs. j23| and |22I is not required for that purpose. Moreover, a possible small 
underestimation of ^ in the broken phase can be controlled by comparing the critical 
amplitude ratio C,^ /C,^, which we calculate below, to the known results for the ratio 
Uf . 

<,exp 



18 




6.1 The Critical Amplitudes of the Correlation Length 



We have determined the correlation length on the critical isotherm for if > on 
lattices with L = 72, 96 and 120. The results are shown in Fig.|Hl We have performed 
fits to these data with the ansatz 



e = eH-'^[i+eiH-''^+e2H] . (75) 

As in the case of the magnetization the first correction-to-scaling term turned out 
to be compatible with zero. Subsequently, we have therefore made fits with = 
in the interval 0.0002 < H < 0.0035. Irrespective of whether we use all data in the 
interval, or only those from the largest lattices or subsets of them we find 

^ = 0.3048(9)fr~''^[l + 9.8(2. 5)fr] , (76) 

with x^/Nf varying between 0.8 and 1.1. The corresponding critical amplitude is 
then 

r = 0.3048(9) . (77) 

The correlation length data for H = below and above the critical temperature are 
shown in Figs. El (a) and (b). They are slighly more fluctuating than the data in Fig. 
IHI Nevertheless, the corresponding critical amplitudes may be determined quite well. 
Here, we have tested again that a fit ansatz without the first correction-to-scaling 
term is sufficient 

^ = ^^\tr[i+^tt] . (78) 
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Figure 9: The correlation length ,^ at if = for L = 48, 72 and 96, (a) below Tc 
versus —t and (b) above Tc versus t. The solid lines show the fits from Eqs. (ffn|l 
and (jHTjl . respectively. 

As it turns out, the best fits in the broken phase (with /Nf ^ 3) are those using 
all points and setting .^g' = 0. We obtain 

^ = 0.2145(15) ■ (-^)"^ that is C = 0.2145(15) . (79) 

The error of ^~ includes also all other fit results, using only the data from the largest 
lattices and/or with a correction term proportional to t. In the symmetric phase the 
correction term is only observable, when the larger t-points are included in the fits. 
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We have therefore first performed fits in the interval 0.002 < t < 0.009 omitting the 
correction term and find for the critical amplitude the value 

^+ = 0.415(1) . (80) 

The subsequent fit with all data leads to 

^ = 0.415(l)r''[l + 0.2(l)t] , (81) 

with again /Nf 3. From the critical amplitudes for if = we derive the ratio 
^+/^~ = 1.935(14). As discussed already, this number is to be compared to the 
results 1.891(22) from [2H1, |22l El and 1.896(10) from EZj for this ratio. We 
deduce from this check, that our ^-values in the broken phase are about 2% lower 
than the true exponential correlation length. The critical amplitudes and 
which we have also determined, should however be unaffected by possible higher 
states [22] • That allows us to calculate the universal ratios Qc and Q2 from Eqs. 
(Uni) and (1201). We obtain 

= 0.326(3) and Q2 = 1.201(10) . (82) 

Our finding for is well in accord with two other results: 0.328(5) from and 
0.324(6) from |2E], it is slightly below the value 0.3315(10) from 01211. The number 
for Q2 is the first one from MC simulations and compares well with 1.195(10) from 
0I2Z1, the result 1.17(2) from [1311221 is somewhat lower. 



7 The Scaling Function of the Correlation Length 

In Fig. ^1 we show the data for the normalized scaling function g^{z) = C.h'^" / g^{0)- 
The determination of the correlation length was sometimes difficult, in particular in 
the region below Tc, because the correlation function data were not always accurate 
enough to identify a reasonable plateau in ^e//- In niost of these cases we could resort 
to the results of smaller lattices, after convincing ourselves from data adjacent in 
H, that no finite size effects were present. The data shown in Fig. ^1 correspond 
to the ones presented in Fig. [HI for the scaling function of the susceptibility, apart 
from a point at z = —10.4. The normalization was calculated from Eq. (jH^ and 
is g^{0) = 0.331(1). As in the case of the magnetization and the susceptibility 
we observe scaling of the correlation length data in a large z-range and an early 
approach to the asymptotic forms calculated from Eq. ()40|) . The shape of g^{z) is 
very similar to that of /^(-z) - both functions have a peak at about the same 2;- value, 
slightly below 2. It is therefore natural to consider the ratio g'^/fx appearing in Eq. 

Taking the square of g^ instead of the function itself leads to a very slowly 
varying function of z. That can be seen from the asymptotic behaviour 

9l/fxi^).=^o.'^(^'y"' ■ (83) 
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Figure 10: The scaling function of the correlation length: g(^{z) = ^h"^" /g^^iO). The 
dashed lines are the asymptotic forms calculated from Eq. ()4()j) . The numbers refer 
to the different J = 1/T- values of the data. 
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Figure 11: The ratio g'^/f^{z) of scaling functions. The dashed lines are the asymp- 
totic forms calculated from Eqs. ()8H|1 and (jH3)- The numbers refer to the different 
J = 1/T-values of the data, the star shows the normalization point ^|//x(0) = S. 
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Here, rju = 0.0211 and 



1 l/(2-r;) 



and 



Uo 



For 77 — »• one obtains an even simpler picture: the ratio gf/ fx attains the constant 
values for (±^) — >■ 00, where r+ = S/Q2 = 4.001 and r^/r~ is very close to 3/4 
(when one takes the "true" value = 1.89), at 2; = one has g^/ fx = ^ = 4.805. In 
Fig.^Jwe show the behaviour of the ratio from the actual data and we compare it to 
the asymptotic forms (jH^ derived from the amplitudes which we have determined. 
As expected, the ratio is slowly varying, apart from the close neighbourhood of 
the critical line 2; = 0. In the hot phase {z > 0) we obviously have a consistent 
behaviour of the data and a clear, early approach to the asymptotic curve, whereas 
in the cold phase {z < 0) the data are more noisy, probably because of the difficulties 
in calculating the correlation length with sufficient precision there. We refrain here 
from a direct comparison with the parametric representation of Fisher et al. jH] for 
^^/2x because we use a slightly different set of critical exponents (those of Ref. 
jni), which is crucial here. Moreover, since we consider an improved instead of the 
standard Ising model the relative normalizations of the variables are unclear. 



8 Summary and Conclusions 



In this paper we have calculated the equation of state for the universality class of the 
three-dimensional Ising model directly from Monte Carlo simulations. In addition 
we have determined the universal scaling functions of the susceptibility and of the 
correlation length. There are essentially three difficulties one encounters in such 
an enterprise: finite size effects, corrections to scaling and the unknown size of the 
critical region, where scaling works and the equation of state describes the critical 
behaviour of the magnetization. In order to minimize the corrections to scaling we 
have studied the improved Ising model as proposed in Ref. We used cubic 

lattices of linear extension L in the range 24 to 120 to prevent finite size effects and 
performed a large number of simulations with non-zero magnetic field H at different 
couplings J = 1/T to identify the critical region. Further simulations at if = 
and at T = Tc served to determine the critical amplitudes of the magnetization, 
the susceptibility and the correlation length. In particular, we have calculated the 
amplitudes B'^ of the magnetization and ^'^ of the correlation length for the first 
time from MC simulations. The critical behaviours of M, x ^ these lines 
in the H, T-plane could indeed be described well without the leading correction-to- 
scaling terms. From the amplitudes, which are different from those of the standard 
Ising model, we have derived the universal amplitude ratios C~^/C~ = 4.756(28), 
Rx = 1.723(13), Qe = 0.326(3) and Q2 = 1.201(10), in accord with other MC and 
analytic results. Only our estimate for = 1.935(14) differs from the corresponding 
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ratio for the exponential correlation length, probably because the .^-values in the 
broken phase were underestimated slightly due to the contribution of higher states 
to the correlation functions. 

As it turned out, the critical region extends only to rather small values in the 
external field if, especially in the region below the critical temperature. We have 
compared our scaling magnetization data to the parametric representation of the 
equation of state by Guida and Zinn- Justin and find excellent agreement 

between theory and numerical experiment. The comparison of the susceptibility data 
to the corresponding scaling function shows a marginal difference in the symmetric 
phase, which however can be explained by the slightly different value for Ry. used 
in the parametrization. The shape of the correlation-length-scaling function is very 
similar to the one of the susceptibility, the peak positions of the two scaling functions 
are coinciding within the error bars. The ratio is therefore a slowly varying 

function without a peak, as expected from earlier parametrizations. 
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